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Abstract — This paper explains how a popular, commercially- 
available software package for solving partial-differential- 
equations (PDEs), as based on the finite-element method (FEM), 
can be configured to calculate, efficiently, the frequencies and 
fields of the whispering-gallery (WG) modes of axisymmetric 
dielectric resonators. The approach is traceable; it exploits 
the PDE-solver's ability to accept the definition of solutions 
to Maxwell's equations in so-called 'weak form\ Associated 
expressions and methods for estimating a WG mode's volume, 
filling factor(s) and, in the case of closed(open) resonators, its 
wall(radiation) loss, are provided. As no transverse approxima- 
tion is imposed, the approach remains accurate even for quasi- 
transverse-magnetic/electric modes of low, finite azimuthal mode 
order. The approach's generality and utility are demonstrated by 
modeling several non- trivial structures: (i) two different optical 
microcavities [one toroidal made of silica, the other an AlGaAs 
microdisk]; (ii) a 3rd-order sapphireiair Bragg cavity; (iii) two 
different cryogenic sapphire WG-mode resonators; both (ii) and 
(iii) operate in the microwave X-band. By fitting one of (iii) to 
a set of measured resonance frequencies, the dielectric constants 
of sapphire at liquid-helium temperature have been estimated. 



I. Introduction 

NON-TRIVIAL electromagnetic structures can be mod- 
eled through computer-aided design (CAD) tools in 
conjunction with programs for numerically solving Maxwell's 
equations. Though alternatives abound [1], [2], [3], the latter 
often use the finite-element method (FEM) [4], [5]. Within 
such a scheme, a problem frequently encountered when at- 
tempting to determine the values of electromagnetic parame- 
ters from experimental data is a lack of traceability: significant 
dependencies between the data, the model's configurational 
settings, and the inferred values of parameters cannot be 
adequately isolated, understood, or quantified. Traceability de- 
mands that both the model's definition and its solution remain 
amenable to complete, explicit description. And, furthermore, 
convenience requires that the representations adopted for this 
purpose be concise -yet wholly unambiguous. 

A. Whispering -gallery -mode resonators 

Certain compact electromagnetic structures support closed 
whispering-gallery (WG) modes. Though elliptical [6] or even 
non-planar ('crinkled' [7] or 'spooled' [8]) WG morpholo- 
gies exhibit advantageous features with respect to certain 
applications, this paper considers only (closed, planar) WG 
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modes with circular trajectories, as supported by axisymmetric 
resonators within the class depicted in Fig. 1, upon which 
various recent scientific innovations [9], [10], [11], [12] are 
based. Many of the current commercial software packages for 




Fig. 1. Generic axisymmetric resonator in cross-section (medial half -plane). 
A dielectric volume (in 3D) or 'domain' (in 2D) is enclosed by an electric 
wall (El) -or one subject to some different boundary condition, as per 
subsubsection III-E.2. This domain comprises several subdomains (Dl, D2, 
and D3), each containing a spatially uniform dielectric (that could be just 
free space). Some of these subdomains (D2 and D3) are bounded internally 
by electric walls (E2, E2^ and E3). The resonator has (optionally) a mirror 
symmetry through its (horizontal) equatorial plane (dashed line Ml); on 
imposing an electric or magnetic wall over this plane, only either the upper 
or lower half of the resonator need be simulated. 



modeling electromagnetic resonators suffer, however, from a 
'blind spot' when it comes to calculating, efficiently (hence 
accurately), such resonators' whispering-gallery modes. The 
popular MAFIA/CST package [13] is a case in point: as Basu 
et al [14] and no doubt others have experienced, it cannot 
be configured to take advantage of a circular WG mode's 
known azimuthal dependence, viz. exp(±iM(/)), where M (an 
integer > 0) is the mode's azimuthal mode order, and (/) the 
azimuthal coordinate. Though frequencies and field patterns 
can be obtained (at least for WG modes of low azimuthal 
mode order), the computationally advantageous reduction of 
the problem from 3D to 2D that the resonator's rotational 
symmetry affords is, consequently, precluded^ ; and the ability 
to simulate high-order WG modes with sufficient accuracy (for 
metrological purposes) is, exasperatingly, lost. 

^ About the best one can do is simulate a 'wedge' [over an azimuthal domain 
A(/) = 7r/(2M) wide] between radial electric and magnetic walls. 
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B. Brief, selected history of WG-mode simulation 

The method of 'separating the variables' provides analytical 
expressions for the WG modes of right-cylindrical uniform 
dielectric cavities (or shells) [15], [16], [17]. By matching 
expressions across certain boundaries, approximate WG-mode 
solutions for composite cylindrical cavities can be obtained 
[18], [19], [9], whose discrete/integer indices (related to 
symmetries) provide a nomenclature [20] for classifying the 
lower-order WG modes of all similar structures. Extensions 
of the basic mode-matching method encompassing spatially 
non-uniform field polarizations have been developed [3]. 

The accurate solution of arbitrarily shaped axisymmetrical 
dielectric resonators requires numerical methods. Apart from 
the finite-element method (FEM) itself, the most developed 
and (thus) immediately exploitable alternatives include (given 
here only for reference -not considered in any greater de- 
tail): (i) the Ritz-Rayleigh variational or 'moment' meth- 
ods [21], [22], [23], (ii) the finite difference time domain 
method (FDTD) [1], [24], and (iii) the boundary-integral 
[2] or boundary-element methods (BEM, including FEM- 
based hybridizations thereof [25]). Zienkiewicz and Taylor 
[4], particularly their table 3.2, indicate various commonalities 
between them (and FEM).^ 

The application of the finite-element method to the solving 
of Maxwell's equations has a history [26], and is now an 
industry [13], [27], [28]; ref. [4] supplies FEM's theoretical 
underpinnings. Though the method can solve for all three of 
a WG mode's field components, the statement of Maxwell's 
(coupled partial-differential) equations in component form can 
be onerous, if not excluded outright by the equation- solving 
software's lack of configurability. With circular whispering- 
gallery modes, the configurational effort can be significantly 
reduced by invoking a so-called 'transverse' approximation 
[14], [29], wherein only a single (scalar) partial-differential 
equation is solved (in 2D). Here, either the magnetic or 
electric field is assumed to lie everywhere parallel to the 
resonator's axis of rotational symmetry (see figure B.l of 
ref. [29]). This approximation is, however, uncontrolled.^ This 
paper demonstrates that, through only a modicum of extra 
configurational effort, the transverse approximation and its 
associated doubts can be wholly obviated. 

A problem that besets the direct application of FEM to the 
solving of Maxwell's equations is the generation of (many) 
spurious solutions [30], [31], associated with the local gauge 
invariance, or 'null space' [31], that is a feature of the 
equations' 'curl' operators. At least two research groups have 
nevertheless developed in-house software tools for calculating 
the WG modes of axisymmetric dielectric resonators, that: 
(i) solve for all field components (i.e. no transverse approxi- 
mation is invoked), (ii) are 2D (and thus numerically efficient) 

^It is remarked parenthetically here that FDTD may be regarded as a 
variant of FEM employing local, discontinuous shape functions. It is perhaps 
also worth acknowledging that, for resonators comprising just a few, large 
domains of uniform dielectric, the boundary-integral methods (based on 
Green functions), which -in a nutshell- exploit such uniformity to reduce 
the problem's dimensionality by one, will generally be more computationally 
efficient than FEM. 

^It is noted paranthetically that basic mode matching [9] also invokes the 
same transverse approximation and is thus equally uncontrolled. 



and (iii) effectively suppress spurious solutions (without detri- 
mental side-effects) [30], [21], [22], [32], [33]. The method 
described in this paper sports these same three attributes. With 
regard to (iii), the approach adopted by Auborg et al was 
to use different finite elements (v/z. a mixture of 'Nedelec' 
and 'Lagrange' -both 2nd order) for different components 
of the electric and magnetic fields; Osegueda et al [32], on 
the other hand, used a so-called 'penalty term' to suppress 
(spurious) divergence of the magnetic field. Stripping away its 
motivating remarks, applications and illustrations, this paper, 
in essence, translates the latter approach into explicit 'weak- 
form' expressions, that can be directly and openly ported to 
any partial-differential equation (PDF) solver (most notably 
COMSOL/FEMLAB [27]) capable of accepting such. 

II. Method of solution 

A. Weak forms 

Scope: The resonators treated below comprise volumes of 
dielectric space bounded by either electric or magnetic walls 
(or both) -see again Fig. 1; the restriction to axisymmetric 
resonators is only invoked at the start of subsection II-B. The 
resonator's dielectric space comprises voids {i.e. free space) 
and pieces of (sufficiently low-loss) dielectric material; its 
(default) enclosing surfaces will generally be metallic, corre- 
sponding to electric walls. When modeling resonators whose 
forms exhibit reflection symmetries, where the magnetic and 
electric fields of their solutions transform either symmetrically 
or antisymmetrically through each mirror plane, perfect mag- 
netic and electric walls can be alternatively imposed over these 
planes to solve for different 'sectors' of solutions. 

The electromagnetic fields within the dielectric volumes of 
the resonator obey Maxwell's equations [34], [17], as they are 
applied to continuous, macroscopic media [35]. Assuming the 
resonator's constituent dielectric elements have negligible (or 
at least the same) magnetic susceptibility, the magnetic field 
strength H will be continuous across interfaces.^ This property 
makes it advantageous to solve for H (or, equivalently, the 
magnetic flux density B = /iH -with a constant global 
magnetic permeability /i), as opposed to the electric field 
strength E (or displacement D). Upon substituting one of 
Maxwell's curl equations into the another, the problem reduces 
to that of solving a (modified) vector Helmholtz equation: 

V X (e-^ V X H) - aV{V • H) + c'^d'^YS./ dt^ = 0, (1) 

subject to appropriate boundary conditions (read on). Here, c 
is the speed of light and the inverse relative permittivity 
tensor; one assumes that the resonator's dielectric elements are 
linear, such that e"^ is a (tensorial) constant -i.e. independent 
of field strength. Providing no magnetic monopoles lurk inside 
the resonator. Maxwell's equations demand that H be free of 
divergence, i.e. V-H = 0. The middle, so-called 'penalty' term 
on the left-hand side of equation 1 acts to suppress spurious 
solutions, for which (in general) V • H ^ 0; it has exactly 

^The method described in this paper could be to extended to treat resonators 
containing dielectrics with different magnetic susceptibihties by setting up 
(within the PDE-solver -i.e. COMSOL) 'coupling variables' at interfaces. 
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the same form as that used by Osegueda et al [32].^ The 
constant ol controls the penalty term's weight with respect 
to its Maxwellian neighbors; a = 1 was taken for every 
simulation presented in sections IV and V. 

Reference [17] (particular section 1.3 thereof) supplies a 
primer on the electromagnetic boundary conditions stated 
here. The magnetic flux density at any point on a (perfect) 
electric wall satisfies B • n = 0, where n denotes the wall's 
surface normal vector. Providing the magnetic susceptibility 
of the dielectric medium bounded by the electric wall is not 
anisotropic, this condition is equivalent to 

H . n = 0. (2) 

The electric field strength at the electric wall obeys 

E X n = 0. (3) 

These two equations ensure that the magnetic(electric) field 
strength is oriented tangential(normal) to the electric wall. 
As is pointed out in reference [32], equation 3 is a so-called 
'natural' (or, synonymously, a 'naturally satisfied') boundary 
condition within the finite-element method -see ref. [4]. 

The boundary conditions corresponding to a perfect mag- 
netic wall (dual to the those for an electric wall) are 

D • n = 0, (4) 

and 

H X n = 0; (5) 

these two equations ensure that the electric displace- 
ment(magnetic field strength) is oriented tangential(normal) 
to the magnetic wall. Again, the latter equation is naturally 
satisfied. 

One now invokes Galerkin's method of weighted residuals 
[4]; ref. [31] provides an analogous treatment when solving 
for the electric field strength (E). Both sides of equation 1 
are multiplied (scalar-product contraction) by the complex 
conjugate of a 'test' magnetic field strength H , then integrated 
over the resonator's dielectric volume. Upon expanding the 
permittivity-modified 'curl of a curl' operator (to extract a 
similarly modified Laplacian operator), then integrating by 
parts (spatially), then disposing of surface terms through the 
electric- or magnetic-wall boundary conditions stated above, 
one arrives (equivalent to equation (2) of reference [32]) at 

/ [(V X H*) - (V X H) - 
Jv ^ 

ce(V • H*)(V • H) + c-^H* • a^H/Qt^] dV = 0, (6) 

where 'Jy' denotes a volume integral over the resonator and 
(V X H*)-(V X H) = E- i=i[^ X H*]ie-.i[V x H],-, where 
e-^^ are the components of the inverse relative permittivity 
tensor. The three terms appearing in the integrand correspond 
directly to the three weak-form terms required to define an 
appropriate model within a partial-differential-equation solver. 

^Though COMSOL can implement mixed ('Nedelec' plus 'Lagrange') finite 
elements [21], it was found that equation I's penalty term (with its weighting 
factor somewhere in the range 0.01 < ol < 10) could, in conjunction with 
2nd-order Lagrange finite elements (applied to all three components of H), 
always satisfactorily suppress the spurious modes. 



Assuming that the physical dimensions and electromagnetic 
properties of the resonator's components are temporally in- 
variant (or at least 'quasi- static'), solutions or 'modes' take 
the form H(r; t) = H(r)exp(— i27r/t), where r is the vector 
of spatial position, t the time, and / the mode's resonance 
frequency. The last, 'temporal' term in equation 6's integrand 
can thereupon be re-expressed as — (c/)^H(r)* • H(r), where 
c = 27r/c; this re-expression reveals the integrand's complete 
dual symmetry between H and H. 

B. Axisymmetric resonators 

The analysis is now restricted to axisymmetric resonators, 
where a system of cylindrical coordinates (see top right Fig. 1), 
aligned with respect to the resonator's axis of rotational sym- 
metry, has components {r, 0, = {'rad(ial)', 'azi(muthal)', 
'axi(ial)'}. The aim is to calculate the resonance frequencies 
and field patterns of the resonator's circular WG modes, whose 
phase varies as exp(iM(/)), with M = {0, 1, 2, ...} the mode's 
azimuthal mode order. ^ Viewed as a three-component vector 
field over (for the moment) a three-dimensional space, the 
time-independent part of the magnetic field strength now takes 
the form 

H(r) = e^^^ { iy^(r, z),i H^{r, z), H\r, z) } (7) 

where an 'i' (= a/(— 1)) has been inserted into the field's az- 
imuthal component to allow, in subsequent solutions, all three 
component amplitudes iJ"^, to be each expressible 
as a real amplitude multiplied by a common complex phase 
factor. The relative permittivity tensor of an axisymmetric 
dielectric material is diagonal with entries (running down the 
diagonal) €diag. = {e_L, e_L, ey}, where e||(e_L) is the material's 
relative permittivity in the axial direction (in the transverse or 
'perpendicular' plane -as spanned by its radial and azimuthal 
directions). 

One now substitutes equation 7 into equation 6's integrand; 
textbooks provide the required explicit expressions for the 
vector differential operators in cylindrical coordinates [16], 
[17]. A radial factor, r, is included here from the volume 
integral's measure: dV = 2TTrdrd(j)dz (the common factor 
of 27r is dropped from all expressions below.) The first, 
'Laplacian' weak term is given by 

(V X H*)- (V X H) = + 5 + ^c)/(6^6||), (8) 

where 

A={ e_L[^^i?^ - M{H'^H'^ + i7^^^) + M^H'^H'^)] 

^e\\M^H'H'}, (9) 

- e\\M{H'Hf + H'Hf), (10) 

^The method does not require M to be large; even modes that are 
themselves axisymmetric, corresponding to M = 0, such as the one shown 
in Fig. 6(b), can be calculated. 
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C = {e^HfHf 



Note that the transformation {n^ 



Ur}, connects 



H:){H^-H:)+HtHt]}. (11) 



Notation: denotes the partial derivative of (the az- 
imuthal component of the magnetic field strength) with respect 
to r (the radial component of displacement), etc.. Here, the 
individual factors and terms have been ordered and grouped 
so as to display the dual symmetry. Similarly, the weak penalty 
term is given by 



a{V ' H*)(V • H) = a{— + E + r F}, (12) 



where 



D = H^'H'' - M{H'^H'' + H'^H'') + M'^H'^H'^, (13) 

E = {H; + - MH"^) 

+ {H^- - MH'^){Hl + HI), (14) 

F^{H:^HI){HI^HI). (15) 

And the temporal weak-form ('dweak') term is given by 

= -c^f r {H'^H'^ + J^^i^^ + H'H'), (16) 

where H^^ denotes the double partial derivative of 
w.r.t. time, etc.. Note that none of the terms on the right- 
hand sides of equations 8 through 16 depend on the azimuthal 
coordinate 0; the problem has been reduced from 3D to 2D. 

C. Axisymmetric boundary conditions 

An axisymmetric boundary's unit normal in cylindrical 
components can be expressed as {71^,0,712} -note vanishing 
azimuthal component. The electric-wall boundary conditions, 
in cylindrical components, are as follows: H • n = gives 

H^Ur + W^n^ = 0, (17) 

and E X n = gives both 

Hl-H^ = (18) 

and 

e^{H'^-H'^M^Hfr)nr-e\\{H'M-Hfr)n, =0. (19) 

When the dielectric permittivity of the medium bounded by 
the electric wall is isotropic, the last condition reduces to 

{H"^ - H^'M + H^r)nr - {H^M - Hfr)n^ = 0. (20) 

The magnetic-wall boundary conditions, in cylindrical com- 
ponents, are as follows: D • n = gives 

{H^M - Hfr)nr + {H'f' - H^'M + Hfr)n^ = 0, (21) 

and H X n = gives both 

H^rir - H^'n^ = (22) 

and 

H"^ = 0. (23) 



equations 17 and 22, and equations 20 with 21. The above 
weak-form expressions and boundary conditions, v/z. equa- 
tions 8 through 23 are the key results of this paper: once 
inserted into a PDE- solver, the WG modes of axisymmetric 
dielectric resonators can be readily calculated. 

III. Postprocessing of solutions 

Having determined, for each mode, its frequency and all 
three components of its magnetic field strength H as a function 
of position, other relevant fields and parameters can be derived 
from this information. 

A. Other fields (related through Maxwell's equations) 

Straightaway, the magnetic flux density B = /iH. As no real 
('^(9^-displacement') current flows within a dielectric, V x 
H(t) = dD{t)/dU thus D = -i(27r/)-^V x H{t). And, E = 
e~^D, where is the (diagonal) inverse permittivity tensor, 
as already discussed in connection with equation 6 above. 

B. Mode volume 

Accepting various caveats (most fundamentally, the problem 
of mode- volume divergence -see footnote 8) as addressed by 
Kippenberg [29], the volume of a mode is defined as [12] 

_JJJ^_, elEpdV 
^"^"^^ ^ max[e|EP] ' ^^"^^ 

where max [...], denotes the maximum value of its functional 
argument, and / / Jh s •••^^ denotes integration over and 
around the mode's 'bright spot' -where its electromagnetic 
field energy is concentrated. 

C. Filling factor 

The resonator's electric filling factor, for a given dielectric 
component (labeled 'diel.'), a given mode, and a given field 
direction, ('dir.' G {radial, azimuthal, axial}), is defined as 



TTidir. 



(25) 



///e|E|2dV ' 

where J J J^.^^ ...dV denotes integration over the component's 
volume and pol. = {_L, ||} for dir. = {radial or azimuthal, 
axial}. The numerators and denominators of equations 24 
and 25 can be readily evaluated using the PDE-solver's post- 
processing features. 

D. Wall loss (closed resonators) 

Real resonators suffer losses that render the Qs of their 
modes finite. The energy stored in a mode's electromag- 
netic field is U = {1/2) J J J /i|HpdV. For axisymmetric 
resonators, the 3D volume integral / / / dV reduces to a 
2D integral / / {27rr)drdz over the resonator's medial half- 
plane. The surface current induced in the resonator's enclosing 
electric wall is (see ref. [17], page 205, for example) = 
nxH; the averaged-over-a-cycle power dissipated in the wall is 
PiQgg = (1/2) / / i^^lHtpdS, where Ht is the tangential (with 



TRACEABLE 2D FINITE-ELEMENT SIMULATION OF WHISPERING-GALLERY MODES 



5 



respect to the wall) component of H, Rg = (7r////cr)^/^ is the 
wall's surface resistivity (see refs. [17], [34]), a the wall's 
(bulk) electrically conductivity,^ and / the mode's frequency. 
For axisymmetric resonators, the 2D surface integral / / dS 
reduces to a ID integral / (27rr)dl around the periphery of the 
resonator's medial (r-z) half-plane. The quality factor, defined 
as 27r/ U/P\oss, due to the wall loss can thus be expressed as 



Owaii = ^A=(47r/Ma)i/% 



(26) 



where A, which has the dimensions of a length, is defined as 
//iHtPdS 



'dl 



(27) 



Again, both integrals (numerator and denominator), hence 
Qwaii itself, can be readily evaluated through the PDE- solver's 
post-processing features. 

E. Radiation loss (open resonators) 

Preliminary remarks: With open whispering-gallery-mode 
resonators (either microwave [38] or optical [10], [12]), the 
otherwise highly localized WG mode spreads throughout free- 
space;^ energy flows away from the mode's bright spot (where 
the electric- and magnetic-field amplitudes are greatest) as 
radiation. The tangential electric and magnetic fields on any 
closed surface surrounding the bright spot constitute, by the 
'Field Equivalence Principle' [39], [40] (as a formalization of 
Huygen's picture), a secondary source of this radiation. 

1 ) Underestimator of loss via retro -reflection: Consider a 
closed, equivalent of the open resonator, with an enclosing 
electric wall in the WG-mode's radiation zone. The wall's 
form is chosen such that -as far as possible- the open 
resonator's radiation propagates (as a predominantly transverse 
and locally plane wave) in a direction that is locally normal 
to the wall. The electric wall then reflects the otherwise open 
resonator's radiation back onto itself -so creating a standing 
wave, i.e. 3. loss-less mode. Through an argument reminiscent 
of Schelkunoff's induction theorem [39], [41], the tangential 
magnetic field of this mode, H^^°^®^, at any point just inside 
the closed resonator's electric wall, can be related to that of the 
corresponding open resonator's radiation, H°^®", at the same 
point, through H°p^" > 2Ht^^°"^^. The radiation loss can be 
evaluated by integrating the cycle-averaged Poynting vector 
over the electric wall, i.e. Prad. = (1/2) / / ^o|H°P^"pdS, 
where zq is the impedance of free space. A bound on the 
mode's radiative Q-factor can thus be expressed as 



Qrad. > (87r//c)A, 



(28) 



^It is here pointed out that, at Hquid-heHum temperatures, the bulk and 
surface resistances of metals can depend greatiy on the levels of (magnetic) 
impurities within them [36], and the text-book /~^/^ dependence of surface 
resistance on frequency is often violated [37]. 

^As understood by Kippenberg [29], this observation implies that the 
support of equation 24's f f _^ ...dV integral (spanning the WG mode's 
bright spot) must be somehow limited, spatially, or otherwise (asymptotically) 
rolled off, lest the integral diverge. [The so-called 'quantization volume' 
associated with the radiation extends to infinity.] 



approaching equality when the WG mode's bright spot lies 
(in effect) at an antinode of the mode's standing wave; A here 



is exactly that given by equation 27, with H = H! 



closed 



only 



cos 



now the enclosing electric wall is in the radiation zone. 

2) Overestimator of loss via outward-going free-space 
impedance matching: A complementary bound can be con- 
structed by replacing the above closed resonator's electric wall 
with one, of the same form, that attempts to match, impedance- 
wise, the open resonator's radiation -and thus absorb it. 
For transverse, locally plane-wave radiation in the radiation 
zone (in free space) sufficiently far from the resonator, the 
required impedance-matching boundary condition on the wall 
is n X H = E — n(E • n),^ where n is the wall's inward- 
pointing normal. Upon differentiating with respect to time and 
using IMaxwell's displacement-current equation, this condition 
can, for a given mode, be generalized to 

(^mix){VxH-n[(VxH)-n]} 

+ sin(6>mix) i cf mode n X H = 0, (29) 

where /mode is the mode's frequency, c = 27r/c as before, 
and ^mix is a 'mixing angle' for impedance matching (with 
respect to an outward-going radiation), 6>mix = 7r/4. Unless 
^mix = A^7r/2 for integer A^, the i = in equation 29 

breaks the hermitian-ness of the matrix that the PDE-solver 
is required to eigensolve, leading to decaying modes with 
complex eigenfrequencies /mode, and corresponding quality 
factors equal [12] to 3?[/mode]/2^[/mode], where and 
^[...] denote real and imaginary parts. Without fine tuning, 
the enclosing wall's shape will not everywhere lie exactly 
orthogonal to the direction of propagation of the WG mode's 
radiation; thus, even for 6>mix = 7r/4, the radiation will not 
be completely absorbed at the wall. A bound on the mode's 
radiative Q-factor can thus be expressed as 



Qrad. ^ ^[/mode]/2^[/mode] 5 



(30) 



approaching equality on perfect absorption (no reflections). 

IV. Example APPLICATIONS 

The source codes and configuration scripts used to imple- 
ment the simulations presented in this and the next section are 
freely available from the author. 



A. 'Sloping -shoulders' cryogenic sapphire microwave res- 
onator [UWA] 

This axisymmetric resonator [42] comprises a piece of 
monocrystalline sapphire mounted (co-axially) within a cylin- 
drical metal can -see Fig. 2(a); the crystal's optical (or 'c') 
axis lies parallel to the resonator's geometric axis. The piece's 

^Note that the direction (polarization) of E or H in the wall's plane is 
not constrained; the two fields need only be orthogonal with their relative 
amplitudes in the ratio of the impedance of free space zq. 

^^The two terms on the left-hand side of equation 29 can be viewed as 
implementing electric- (cf. equation 3) and magnetic-wall (cf. equation 5) 
boundary conditions, respectively. The (composite) boundary condition can be 
continuously adjusted between these two cases by varying the mixing angle 
^mix- [Parenthetically: ^mix = — 7r/4 corresponds to impedance matching 
inward-coming (as opposed to an outward-going) radiation.] 
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TABLE I 

Electric filling factors for the WGEi4,o,o mode of UWA's 
sloping-shoulders resonator 



Fig. 2. UWA's sloping-shoulder cryogenic sapphire resonator: (a) medial 
cross-section through its geometry; the grey(white) shading corresponds to 
sapphire(vacuum); Si and S2 indicate the sapphire piece's upper and lower 
'shoulders', (b) mesh of the resonator's model structure generated by the 
FEM-based PDE-solver; for clarity, only every other meshing line is drawn 
[i.e. (b) displays the 'half-mesh']; within (c) and (d), the (logarithmic) grey 
scale reflects the absolute value of the vectorial magnetic H and electric E 
fields, respectively; white arrows indicate the magnitude and direction of each 
field's medial component. 



sloping shoulders (Si and S2 ibid.) make accurate simulation 
via mode-matching less straightforward. The resonator's form, 
as encoded into the PDE-solver, is taken from figure 3 of 
ref. [9]^\ with the piece's outer diameter, the length of its 
outer axial sidewall, the axial extent of each sloping shoulder, 
and the radius of each of its two spindles equal to, at liquid- 
helium temperature (i.e. including cryogenic shrinkages -see 
section V) 49.97, 19.986, 4.996, and 19.988 mm, respectively. 
The sapphire crystal's cryogenic permittivities were taken to 
be {e^, ey} = {9.2725, 11.3486}, as given in ref. [20]. Since 
the sapphire piece and its surrounding metal can do not share 
a common transverse ('equatorial') mirror plane, the speeding 
up of the simulation through the imposition of a magnetic or 
electric wall on such a plane (so halving the 2D region to be 
analyzed) is precluded. 

Fig. 2(b) displays the FEM-based PDE-solver's meshing of 
the model resonator's structure; in COMSOL's vernacular^^, 
the mesh comprises 7296 base-mesh elements and 88587 
degrees of freedom ('DOE'). It took typically 75 seconds, to 
obtain the resonator's lowest (in frequency) 16 modes, for a 
single, given azimuthal mode order M, at [with respect to 
Eig. 2(b)] full mesh density, on a standard, 2004- vintage per- 

^^It is remarked here that the drawn shape of the sapphire piece in figure 3 
of ref. [9] is not wholly consistent with its given dimensions: its outer axial 
sidewall is too long and the slope of its shoulders too slight. 

^^The size/complexity of a finite-element mesh is quantified, within COM- 
SOL Multiphysics, by (i) the number of elements that go to compose its 
so-called 'base mesh' and (ii) its total number of degrees of freedom (DOF) 
-as associated with its so-called 'extended mesh'. 



i^dir. 


radial 


azimuthal 


axial 


sapphire 


0.80922 


0.16494 


7.016 X 10-^ 


vacuum 


0.01061 


8.0533 X 10-3 


1.6543 X 10-"^ 



sonal computer (2.4 GHz, Intel Xeo CPU), without altering the 
PDE-solver's default eigensolver settings. With the azimuthal 
mode order set at M = 14, the model resonator's WGE14 0,0 
mode was found to lie at 11.925 GHz, to be compared with 
11.932 GHz found experimentally [9]. 

Wall loss: This mode's characteristic length A, as evaluated 
with respect to the can wall's enclosing surface, was deter- 
mined to be 2.6 X 10^ m. Based on ref. [37], one estimates 
the surface resistance of copper at liquid-helium temperature 
to be 7 X 10~^ ft per square at 1 1.9 GHz, leading to a wall-loss 
Q of 3.5 X 10^^ for the WGEi4,o,o mode. 
Filling factor: Using equation 25, the electric filling factors for 
the WGE14 0,0 mode were evaluated. These factors, presented 
in TABLE I above, are in good agreement with those labeled 
'EE' in Table IV of ref. [9], that were obtained via a wholly 
independent EEM simulation of the same resonator. 

B. Toroidal silica optical resonator [ Caltech ] 



a: 











j 


D ^^-^^ 







b: 

Fig. 3. (a) Geometry (medial cross-section) and dimensions of a model 
toroidal silica microcavity resonator -after ref. [11]; the torus' principal 
diameter D = 16 /im and its minor diameter d = 3 iim; the central vertical 
dashed line indicates the resonator's axis of continuous rotational symmetry, 
(b) False-color surface plot of the (logarithmic) electric-field intensity |E|^ 
over the dashed box appearing in (a) for this resonator's TEp=i,rn=93 
whispering-gallery mode; white arrows indicate the electric field's magnitude 
and direction in the medial plane. 

The resonator modeled here, based on ref. [11], comprises 
a silica toroid, supported above a substrate by an integral 
interior 'web'; its geometry is shown in Eig. 3(a). The toroid's 
principal and minor diameters are {D^d} = {16,3} /im, 
respectively. The silica dielectric is presumed to be wholly 
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isotropic (i.e., no significant residual stress) with a relative 
permittivity of egii. = 2.090, corresponding to a refractive 
index of tIsh. = y'esii. = 1.4457 at the resonator's op- 
erating wavelength (around 852 nm) and temperature. The 
FEM model's pseudo-random triangular mesh covered an 8- 
by-8 /im square [shown in dashed outlined on the right of 
Fig. 3(a)], with an enhanced meshing density over the silica 
circle and its immediately surrounding free-space; in total, 
the mesh comprised 5990 (base-mesh) elements, with DOF = 
36279. Temporarily adopting Spillane et aPs terminology, the 
resonator's fundamental TE-polarized 93rd-azimuthal-mode- 
order mode (where 'TE' here implies that the polarization of 
the mode's electric field is predominantly aligned with the 
toroid's principal axis -not transverse to it) was found to have 
a frequency of 3.532667 x 10^^ Hz, corresponding to a free- 
space wavelength of A = 848.629 nm (thus close to 852 nm). 
Using equation 24, this mode's volume was evaluated to be 
34.587 /im^; if, instead, the definition stated in equation 5 
of ref. [11] is used, the volume becomes 72.288 /im^ -i.e. a 
factor of rig -J greater. These two values straddle the volume 
of ~55 fim^, for the same dimensions of silica toroid and 
(TE) mode-polarization, as inferred by eye from figure 4 of 
ref. [11]. 

C. Conical microdisk optical resonator [ Caltech ] 

The mode volume can be reduced by going to smaller 
resonators, which, unless the optical wavelength is commen- 
surately reduced, implies lower azimuthal mode orders. The 



a: 



b: 




Fig. 4. (a) Geometry (medial cross-section) and (half-)meshing of model 
GaAlAs microdisk resonator -after ref. [12]; the disk's median diameter is D 
= 2.12 /im and its thickness (axial height) t = 255 nm; its conical external 
sidewall subtends an angle = 26° to the disk's (vertical) axis; for clarity, 
only every other line of the true (full) mesh is drawn. The modeled domain 
in the medial half -plane is a rectangle stretching from 0.02 to 1.5 /xm in 
the radial direction and from -0.5 to +0.5 /im in the axial direction, (b) 
False-color surface plot of the (logarithmic) electric-field intensity |Ep for 
the resonator's TEp=i,rn=ii mode at A = 1263.6 nm; again, white arrows 
indicate the electric field's magnitude and direction in the medial plane. 

model 'microdisk' resonator analyzed here, as depicted in 
Fig. 4(a), duplicates the structure defined in figure 1(a) of 



Srinivasan et al [12]; the disk's constituent dielectric (alternate 
layers of GaAs and GaAlAs) is approximated as a spatially 
uniform, isotropic dielectric, with a refractive index equal 
to n = 3.36. The FEM-modeled domain comprised 4928 
quadrilateral base-mesh elements, with DOF = 60003. Adopt- 
ing the same authors' notation, the resonator's TEp=i^rn=ii 
whispering-gallery mode, as shown in Fig. 4(b), was found 
at 2.372517 x 10^"^ Hz, equating to A = 1263.6 nm; for 
comparison, Srinivasan et al found an equivalent mode at 
1265.41 nm [as depicted in their figure 1(b)]. It is pointed 
out here that the white electric-field arrows in Fig. 4(b) [and 
also, though to a lesser extent, in Fig. 3(b)] are not perfectly 
vertical -as the transverse approximation taken in references 
[9], [29] would assume them to be; the quasi-ness of the true 
mode's transverse-electric polarization is thus revealed. 
Mode volume: Using equation 24, but with the mode excited 
as a standing-wave (doubling the numerator while quadrupling 
the denominator), the mode volume is determined to be 
0.1484 X fim^ 0^ 2.79(A/n)^, still in good agreement with 
Srinivasan et aPs ~2.8(A/n)^. 

Radiation loss: The TEp^i^^^n mode's radiation loss was es- 
timated by implementing both the upper- and lower-bounding 
estimators described in subsection III-E. Here, the microdisk 
and its mode were modeled over an approximate sphere, 
equating to a half-disk in 2D (medial plane). The half- 
disk's diameter was 12 /im and different electromagnetic 
conditions were imposed on its semicircular boundary-see 
Fig. 5.^^ With an electric-wall condition (i.e. equations 2 and 

^ m 

■F p7 V 

a: ^^^^ b: ^^^^ c: ^^^^ 

Fig. 5. Radiation associated with the same [TEp=i,rn=ii, A = 1263.6 nm] 
whispering-gallery mode as presented in Fig. 4 (all three maps use the same 
absolute false-color scale): (a) standing- wave (equal outward- and inward- 
going) radiation with the outer semicircular boundary set as a magnetic wall; 
(b) the same but now with the boundary set as an electric wall; (c) somewhat 
traveling (more outward- than inward-going) radiation with the boundary's 
impedance set to that of an outward-going plane- wave in free space (and with 
the normal magnetic field constrained to vanish). That (c)'s radiation field is 
somewhat dimmer than (b)'s is consistent with the two different estimates of 
the resonator's radiative Q corresponding to (b) and (c) [see text]. 

3 or, equivalently, 17 and {18, 20}) imposed on the half- 
disk's entire boundary [as per Fig. 5(b)], the right-hand of 
equation 28 was evaluated. And, with the E x n = condi- 
tion (viz. equation 3) on the boundary's semicircle replaced 
by the outward-going-plane- wave(-in-free- space) impedance- 
matching condition (viz. equation 29 with 6>mix = 7r/4), while 

^^It is acknowledged that, in reality, the microdisk's substrate would occupy 
a considerable part of half-disk's lower quadrant. 
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the H • n = condition (equation 2) is maintained, the right- 
hand side of equation 30 was evaluated for the radiation pattern 
displayed in Fig. 5(c). For a pseudo-random triangulation mesh 
comprising 4104 elements, with DOF = 24927, the PDF 
solver took, on the author's office computer, 6.55 and 13.05 
seconds, corresponding to Figs. 5(b) and (c), respectively^"^, to 
calculate 10 eigenmodes around 2.373 x 10^^ Hz, of which the 
TFp^i^^^ii mode was one. Together, the resultant estimate 
on the TFp=i,rn=ii mode's radiative-loss quality factor is 
(1.31 < Qrad. < 3.82) X 10^, to be compared with the estimate 
of 9.8 X 10^ (at 1265 nm) reported in table 1 of ref. [12]. 



D. Distributed-Bragg-reflector microwave resonator 

The method's ability to simulate axisymmetric resonators 
of arbitrary cross-sectional complexity is demonstrated here 
by simulating the 10-GHz TFqi mode of a distributed-Bragg- 
reflector (DBR) resonator as analyzed by Flory and Taber 
(F&T) [43] through mode matching. The resonator's model 




a: 

Fig. 6. (a) Geometry (medial cross-section) of a 3rd-order distributed- 
Bragg-reflector resonator within a cylindrical metallic can (hence electric 
interior walls -represented by a solid black rectangle); as per ref. [43], the 
can's interior diameter is 10.98 cm and its interior height is 13.53 cm; the 
horizontal and vertical grey (or pink -in color reproduction) stripes denote 
cylindrical plates and barrels of sapphire; white rectangles correspond to right 
cy Under s/annuli of free-space; the vertical arrow indicates the resonator's axis 
of rotational symmetry, with which the sapphire crystal's c-axis is aligned; 
a magnetic wall is imposed over the resonator's equatorial plane of mirror 
symmetry (dashed horizontal line; cf. Ml in Fig. 1). (b) False-color plot of the 
(logarithmic) electric-field intensity |Ep over the top-right medial quadrant 
of the rotationally invariant (M = 0) TEqi mode; note how the mode is 
strongly localized within the resonator's central cylinder of free-space 

geometry was generated through an auxiliary script written 
in MATLAB. Its corresponding mesh comprised 5476 base- 
mesh elements, with 66603 degrees of freedom (DOF), with 
8 edge vertices for each ~A/4 layer of sapphire. Based on 
ref. [44] 's quartic fitting polynomials, the sapphire crystal's 
dielectric permittivities, at a temperature T = 300 K, were set 
to ex = 9.394 (consistent with ref. [43]) and ey = 11.593. 
The TEoi mode shown in Fig. 6(b) was found to lie at 
10.00183 GHz, in good agreement with F&T's 'precisely 
10.00 GHz'. Using equation 25, the mode's electric filling 
factor for the resonator's sapphire parts was 0.1270, which, 

^^The complex arithmetic associated with the impedance-matching bound- 
ary condition meant that the PDF solver's eigen-solution took approximately 
twice as long to run with this condition imposed -as compared to pure electric- 
or magnetic- wall boundary conditions that do not involve complex arithmetic. 



assuming an (isotropic) loss tangent of 5.9 x 10~^ as per 
F&T, corresponds to a dielectric-loss Q factor of 1.334 mil- 
lion. Through equations 26 and 27, and assuming a surface 
resistance of 0.026 Vt per square as per F&T, the wall-loss 
Q factor was determined to be 29.736 million, leading to 
a composite Q (for dielectric and wall losses operating in 
tandem) of 1.278 million. These three Q values are consistent 
with F&T's stated (composite) Q of '1.3 million, and limited 
entirely by dielectric losses'. 



a: 




b: 



c: 



Fig. 7. (a) Close-up of NPL's cryo-sapphire resonator, with the main 
body of its outer copper can removed. The resonator's chamfered HFMFX 
sapphire ring has an outer diameter of ~46.0 mm and an axial height of 
~25.1 mm. The ring's integral interior 'web', 3mm thick, lies parallel to, 
and is centered (axially) on, the ring's equatorial plane, and is supported 
through a central copper post. Optical refraction at the ring's outer surface 
falsely exaggerates its internal diameter. Above the ring are two loop probes 
for coupling electromagnetically to the resonator's operational whispering- 
gallery mode, (b) geometry of the resonator in medial cross-section; pink/grey 
indicates sapphire, white free space; bounding these dielectric domains, and 
shown as solid black lines, are copper surfaces belonging to the resonator's can 
and web-supporting post [the dashed vertical line (lower right) runs along the 
resonator's cylindrical axis (r = 0)]. (c) false-color map (logarithmic scale) 
of |Hp for the resonator's llth-azimuthal-mode-order fundamental quasi- 
transverse-magnetic (Nlii in ref. [7]'s notation) whispering-gallery mode at 
9.146177 GHz (simulated), as detailed on the 6th row of TABLF II. The 
white arrows indicate the magnitude and direction of this mode's electric 
field strength (E) in the medial plane. 



V. Determination of the permittivities of 

CRYOGENIC SAPPHIRE 

The method is here applied to determine the two dielectric 
constants of monocrystalline sapphire (HEMEX grade [45]) 
at 4.2 K from a set of experimental data, listed in the four 
right-most columns of TABLE II, and corresponding to the 
resonator whose innards are shown in Fig. 7(a). Allowance 
was made for the shrinkages of the resonator's constituent 
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TABLE II 

NPL'S CRYOGENIC SAPPHIRE RESONATOR: SIMULATED AND 
EXPERIMENTAL WG MODES COMPARED 
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^the nomenclature of ref. [7] is used for this column. 
^fuU width half maximum (-3 dB). 

^the difference in frequency between the orthogonal pair of standing-wave 
resonances (akin to a 'Kramers doublet' in atomic physics) associated with the 
WG mode; the experimental parameters stated in other columns correspond 
to the strongest resonance (greatest S'21 at line center) of the pair. 



materials from room to liquid-helium temperature^^ and the 
values of sapphire's two dielectric constants (e^ and ey) were 
initially set equal to those specified in ref. [20]. Fig. 7(b)'s 
geometry was meshed with quadrilaterals over the medial half- 
plane, with 8944 elements in its base mesh, and with DOF = 
108555. For a given azimuthal mode order M, calculating the 
lowest 16 eigenmodes took around 3 minutes on the author's 
office PC (as previously specified). With Fig. 8 as a guide, 
each of the 16 experimental resonances was identified to a 
particular simulated WG mode, as specified in the 4th column 
of TABLE II, lying near to it in frequency; these identifications 
were influenced by requiring that the measured attributes 
(e.g. the FWHM linewidths) of the experimental resonances 
belonging -as per their identifications- to the same 'family' 
of WG modes (e.g. SI, or N2) varied smoothly with M. Filling 
factors were then calculated to quantify the differential change 
in the frequency of each identified mode with respect to 
and ey. The two dielectric constants were then adjusted to 
minimize the variance of the residual (simulated-minus- 
measured) frequency differences. The resultant best-fit values, 
to which the simulated data occupying the three left-most 
columns in TABLE II correspond, were: 



= 9.285 (±0.010); 
en = 11.366 (±0.010). 



(31) 
(32) 



^^By integrating up published linear- thermal-expansion data (viz. Ta- 
ble 4/TABLE 1 of ref. [46]/[47]), sapphire's two cry o- shrinkages were 
estimated to be (1.0 - 7.21 x 10"^) and (1.0 - 5.99 x 10""^) for directions 
parallel and perpendicular to the sapphire's c-axis, respectively. From Table F 
at the back of ref. [48], the cryo-shrinkage of (isotropic) copper was taken to 
be (1.0 - 3.26 x 10-^). 
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Fig. 8. Plot used to aid the identification of experimental with simulated WG 
modes. Solid horizontal lines (16 in total) indicate the center frequencies of 
the former. Solid circles indicate the identification of a simulated mode with an 
experimental one (the difference in their frequencies corresponds to much less 
than a circle's radius in all cases); hollow circles indicate simulated modes that 
were not identified with any experimentally measured one. Quasi-transverse- 
magnetic (q-TM) and quasi-transverse-electric (q-TE) WG modes of the same 
family are joined by (blue-)dashed and (red-)dotted lines respectively; a few 
of the lowest-lying mode families are labeled using standard notation [7]. 



The nominal error assigned to each reflects uncertainties in the 
identifications of certain experimental resonances, each lying 
almost equally close (in both frequency and other attributes) to 
two or more different simulated WG modes. Errors resulting 
from a finite meshing density [33], or those from the finite 
dimensional/geometric tolerances to which the resonator's 
shape was known, were estimated to be small in comparison. 
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